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Abstract 


The classical multi-spectral Atmospheric Correction (AC) algorithm is inadequate for the new generation of 
spaceborne hyperspectral sensors such as NASA’s first hyperspectral Ocean Color Instrument (OCI) onboard 
the anticipated Plankton, Aerosol, Cloud, ocean Ecosystem (PACE) satellite mission. The AC process must 
estimate and remove the atmospheric path radiance contribution due to the Rayleigh scattering by air 
molecules and scattering by aerosols from the measured top-of-atmosphere (TOA) radiance, compensate for 
the absorption by atmospheric gases, and correct for reflection and refraction of the air-sea interface. In 
this work, we present and evaluate an improved AC for hyperspectral sensors developed within NASA’s 
SeaWiFS Data Analysis System software package (SeaDAS). The improvement is based on combining the 
classical AC approach of multi-spectral capabilities to correct for the atmospheric path radiance, extended 
to hyperspectral, with a gas correction algorithm to compensate for absorbing gases in the atmosphere, 
including water vapor. The SeaDAS-hyperspectral version is capable of operationally processing the AC 
of any hyperspectral airborne or spaceborne sensor. The new algorithm development was evaluated and 
assessed using the Hyperspectral Imager for Coastal Ocean (HICO) scenes collected at the Marine Optical 
BuoY (MOBY) site, and other SeaWiFS Bio-optical Archive and Storage System (SeaBASS) and AERosol 
Robotic NETwork - Ocean Color (AERONET-OC) coastal sites. A hyperspectral vicarious calibration was 
applied to HICO, showing the validity and consistency of HICO’s ocean color products. The hyperspectral AC 
capability is currently available in SeaDAS to the scientific community at https://oceancolor.gsfc.nasa.gov/. 
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1. Introduction 


Hyperspectral ocean color remote sensing is at an early stage with respect to algorithm development for 


atmospheric correction and the retrieval of ocean color products (“Pre-Aerosol, Clouds, and ocean Ecosys- 
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tem (PACE) Mission Science Definition Team Report,” 2012). The quantification of global chlorophyll 
concentration [Chl] and other inherent optical properties (IOPs), such as the absorption and backscattering 
coefficients of hydrosols, are derived using multispectral spaceborne radiometers such as NASA’s Moderate 
Resolution Imaging Spectroradiometer (MODIS) (Esaias et al., 1998), currently operating from the Aqua 
spacecraft (MODISA), the Sea-viewing Wide Field-of-view Sensor (SeaWiF'S) that operated from 1997 to 
2010 (McClain et al., 2004) and the Visible Infrared Imaging Radiometer Suite (VIIRS) currently operating 
onboard the Suomi National Polar-orbiting Partnership (Suomi NPP) weather satellite (Murphy et al., 2001). 
Additionally, the European Space Agency (ESA) provided multi-spectral ocean color data from the Medium 
Resolution Imaging Spectrometer (MERIS) operated onboard Envisat, from 2002 to 2012 (Rast et al., 1999). 
Meanwhile, the Ocean and Land Colour Instrument (OLCI) currently operating onboard SENTINEL-3 plat- 
form is the successor to MERIS (Donlon et al., 2012). Utilizing radiometry in a few discrete bands that are 
spectrally located to detect variations in ocean constituents, these sensors have been successful in retrieving 
bulk concentrations of the phytoplankton pigment chlorophyll-a or associated phytoplankton absorption co- 
efficients (Morel and Maritorena, 2001; O’Reilly et al., 1998; Werdell and Bailey, 2005), but the detection of 
specific phytoplankton species or functional groupings using multispectral sensors has been challenging due 
to lack of sufficient spectral resolution to enable distinction (Gordon et al., 2001; Nair et al., 2008). Since 
marine environments are complex ecosystems typically supporting a diverse community of phytoplankton 
species, the link between biodiversity, ecosystem structure, and ecological function is unresolved in current 
ocean color remote sensing science. It was in part to address this knowledge gap that NASA conceived 
of the Plankton, Aerosol, Cloud, ocean Ecosystem (PACE) mission. Onboard the PACE satellite, to be 
launched in 2022, the proposed hyperspectral Ocean Color Instrument (OCI) is being developed to retrieve 
marine reflectance with a continuous spectral sampling of 5-nm resolution from the ultraviolet (UV) to the 
near infrared (NIR) (320 to 900 nm) with six additional short wave infrared (SWIR) channels. Potential 
ocean color products from PACE will include phytoplankton functional types (PFTs), marine net primary 
productivity, and marine carbon stocks, in addition to atmospheric products such as aerosol and cloud prop- 
erties (“Pre-Aerosol, Clouds, and ocean Ecosystem (PACE) Mission Science Definition Team Report,” 2012; 
Werdell, 2016). Although the additional spectral bands of PACE/OCI are expected to significantly advance 
ocean color science, challenges arise for atmospheric correction (Gao et al., 2000, 2009). The AC process 
must estimate and remove the atmospheric path radiance contribution due to the Rayleigh scattering by air 
molecules and scattering by aerosols from the measured top-of-atmosphere (TOA) radiance, compensate for 
the absorption by atmospheric gases, and correct for reflection and refraction of the air-sea interface. The 
residual signal is attributed to the light upwelling from beneath the ocean surface, and is typically <15% 


of the TOA radiance, depending on the spectral range and ocean brightness. This water-leaving signal is 


then used to derive ocean color products assuming that an accurate atmospheric correction was applied 
(Ahmad et al., 2010; Gordon and Wang, 1994a). Multispectral ocean color missions typically avoid light 
detection at the spectral regions that are most affected by absorbing gases, especially water vapor which is 
highly variable spatially and temporally. Current ocean color atmospheric correction algorithms compensate 
mainly for ozone (O3) and nitrogen dioxide (NO2), which absorb light in broad spectral bands from the 
UV and throughout the visible spectral range (Ahmad et al., 2007). Spectral bands in the red and NIR (> 
600 nm) are specifically located to avoid water vapor absorption features that would otherwise complicate 
the derived correction for aerosol scattering and contaminate the derived water-leaving signal (Gordon and 
Wang, 1994a). To make full use of the additional spectral bands of hyperspectral instruments for ocean color, 
the effect of water vapor must be estimated and corrected. The complexity of the atmospheric water vapor 
profile, the spectrally variable nature of the absorption features, and the spatial heterogeneity of the water 
vapor concentration are factors that contribute to this challenge (Gao et al., 2000, 2009). Meanwhile, other 
absorbing gases, discussed in more details in section 2.2, are either typically immutable and homogenous 
or their variability is well quantified based on climatological models or ancillary sources (Dee et al., 2011; 
Derber et al., 1991; Kanamitsu et al., 2002). In this work, we provide a hyperspectral remote sensing applica- 
tion, including AC and ocean color retrieval, based on the Hyperspectral Imager for Coastal Ocean (HICO). 
HICO (Table 1) is a hyperspectral imaging radiometer that operated onboard the international space station 
(ISS) from 2009 to 2014, capturing over 10,000 scenes over the globe (Corson et al., 2008, 2010; Korwan et 
al., 2009a; Lucke et al., 2011). HICO measured light with a spectral coverage from 353 nm to 1080 nm with 
a 5.7 nm spectral resolution. It has a pointing capability in the cross-track direction. At the nadir looking 
direction, the spatial resolution is 90 m. HICO collected one scene per orbit of size 50 x 200 km that was 
scheduled weekly by the science team, with scenes mostly collected over coastal regions to derive products 
such as water clarity, benthic types, and bathymetry. HICO provided adequate radiometric performance to 
support ocean color applications in these coastal regions, where high concentrations of phytoplankton and 
suspended sediments result in strong water-leaving radiance in the visible regime, but scenes collected over 
darker, open ocean regions suffer from the relatively low signal to noise ratio (SNR), especially in the green 
to NIR regime (Korwan et al., 2009a; Lucke et al., 2011). 

The work presented here will provide a detailed description of a hyperspectral atmospheric correction al- 
gorithm that combines the established capabilities of NASA’s standard AC approach mainly for aerosols 
and Rayleigh scattering (L2GEN), currently distributed in SeaDAS (“SeaWiFS Data Analysis System 
(SeaDAS),” 2017), with the absorbing gas corrections found in ATmospheric REMoval code (ATREM), 
discussed in details in section 2.2 (Gao et al., 1993b). This enhanced version of L2GEN extends the multi- 


spectral Rayleigh/aerosol correction heritage approach (Ahmad et al., 2010; Gordon and Wang, 1994a) to the 


Table 1: HICO sensor and data characteristics. 


Platform International Space Station (ISS) 
Operation lifetime 2009-2014 

Orbit repeat time/period 3 days/90 minutes 

Scene size (km) 50 x 200 

Pixel size (m) ~100 

Wavelength (nm) 353 - 1080 (128 bands) 

Spectral resolution (nm) 5.7 

Spectral FWHM (nm) 10 (< 745 nm), 20 (>745 nm) 
Sensor type Offner Spectrometer 
Signal-to-noise ratio (SNR) > 200:1 assuming 5% surface albedo 
Polarization sensitivity <5% 


hyperspectral domain. The novelty of this work is that we have implemented and assessed an operationally 
viable algorithm for hyperspectral ocean color retrieval, that follows the well-established, well validated, 
heritage approach used by NASA for all global OC sensors. Water vapor correction is a significant addition 
to that heritage process, and necessary for meaningful OC retrieval in some spectral regions. The most im- 
portant improvement related to HICO specifically is that HICO’s AC processing has always been dependent 
on in-situ aerosol observations (such as aerosol optical depth, and assumed model). Meanwhile, in this work, 
we provide the community with atmospherically corrected data that is not limited by in-situ observation. 
The AC for HICO is completely automated such that the correction is inferred from HICO measurements 
itself as opposed to in-situ measurements. The general improvement is that this work shows the first op- 
erational hyperspectral AC and water-leaving reflectance retrieval of ocean color sensors. The operational 
algorithm can be the backbone and a starting point for the anticipated PACE mission. We also provide the 
space-borne hyperspectral data for the community to better understand the complication of hyperspectral 
remote sensing of ocean color and develop appropriate algorithms for such capabilities. We also present 
a system-level vicarious calibration to improve the radiometric stability of the HICO sensor and minimize 
bias in the atmospheric correction process (Franz et al., 2007). To assess the performance of the resulting 
water-leaving radiance retrievals, match-ups to in-situ radiometry from SeaBASS and AERONET-OC are 
presented and assessed. Finally, we show some case study results of ocean color retrievals over coastal regions 


in the Chesapeake Bay and the Gulf of Mexico and compare them with co-incident MODISA retrievals. 


2. Approach 


2.1. Atmospheric correction 
The L2GEN (L2 generator), within the SeaDAS software package, is the multi-sensor Level-1 to Level-2 
processing code developed and maintained by NASA’s Ocean Biology Processing Group (OBPG) that is 


capable of retrieving ocean color products from TOA radiances for a host of multispectral sensors. L2GEN 


supports multiple atmospheric correction methods and variations that can be applied to a variety of ocean 
color sensors (Ahmad et al., 2010; Bailey et al., 2010; Gordon and Wang, 1994a; Ruddick et al., 2000; Wang 
et al., 2009). NASA’s heritage atmospheric correction approach assumes the additive property of different 
contributing components in the atmosphere-ocean (AO) system to the top of atmosphere radiance, Lt, which 


can be approximated as: 


Dy (A) = (L(A) + La(A) + Lra(A) + tA) L 5 (A) + TA) Lg (A) + t(A) Lu (A)) X Ty (A) (1) 


where is a sensor spectral band wavelength, L,(A) is the radiance due to multiple scattering by air 
molecules in the absence of aerosols (i.e. Rayleigh scattering), L,(A) is multiple scattering by aerosols in the 
absence of Rayleigh, and L,(A) is the radiance contribution due to Rayleigh-aerosol interactions. T,(A), 
T(A) and t(A) are the gaseous transmittance and the direct and diffuse transmittance, respectively, from Sun 
to surface to the sensor. L,(A) is the surface glint in the specular direction and L(A) is the contribution from 
whitecaps and foam on the surface that is transmitted to the sensor field of view. L(A) is the water-leaving 
radiance that is to be retrieved after the atmospheric correction process, by removing and compensating 
for all other components in Eq. (1). Pre-computed look-up tables (LUT) generated using vector radiative 
transfer (VRT) simulations are used to model the Rayleigh and aerosol contributions (Ahmad et al., 2010; 
Ahmad and Fraser, 1982). The Rayleigh scattering term, which is the dominant contribution over the visible 
regime, is simulated for a wide range of solar and viewing geometries (Ahmad et al., 2010). The hyperspectral 
optical properties used in the VRT simulations, such as the extraterrestrial solar irradiance, Rayleigh optical 
thickness, and the depolarization factors, were optically weighted for the sensor-specific spectral response 
function (Bodhaine et al., 1999; Thuillier et al., 2003). The hyperspectral Rayleigh LUTs are also a function 
of wind speed, to account for secondary AO coupling effects (i.e., the multiple scattered light interactions 
between the atmosphere and the ocean surface), and the inferred Rayleigh scattering coefficient is adjusted 
for surface pressure and terrain height (to support retrievals over inland water masses) (Ahmad and Fraser, 
1982; Bodhaine et al., 1999). Glint and whitecap radiance contributions are modeled based on ancillary wind 
speed and observing geometry (Cox and Munk, 1954; Frouin et al., 1996; Gordon and Wang, 1994b; Koepke, 
1984). The aerosol radiometric contribution is estimated using the method of Gordon and Wang, 1994 
(Gordon and Wang, 1994) with aerosol models developed by Ahmad et al., 2010 (Ahmad et al., 2010), and 
coupled with an iterative bio-optical model to separate the scattering contributions from aerosols and water- 
column constituents in turbid (high-scattering) waters (Bailey et al., 2010). The Gordon and Wang approach 
uses a pair of bands in the NIR or SWIR to derive the spectral slope in measured aerosol reflectance, which 


is then compared to the pre-computed aerosol LUTs representing a range of aerosol models (varying aerosol 


optical properties). The method finds the closest set of aerosol models to the observation, and a weighting 
scheme is then applied to the models that enables the radiance contribution to be extrapolated to the visible 
part of the spectrum. Thus, the aerosol term and the Rayleigh-aerosol interaction contributions can be 
estimated and removed from all spectral bands. The pair of bands in the NIR or SWIR should be sensitive 
to variations in aerosol type and concentration (optical thickness) and avoid contamination from gaseous 
absorption, especially water vapor. Spectral calibration problems in HICO, such as second-order diffraction 
effects on the optical diffraction grating system, contribute to severely degraded radiometric performance 
in the wavelength range > 800 nm, rendering this part of the spectrum unreliable for determination of 
aerosol contribution (Corson et al., 2008; Gao et al., 2012; Gao and Li, 2010; Korwan et al., 2009b). The 
remaining spectral windows in the red/NIR part of the spectrum are few, due to contamination by water 
vapor and oxygen absorption. To minimize these influences, the spectral band pair at 747 nm and 787 nm was 
selected for the aerosol correction. In bright waters, those two bands will be impacted by the non-negligible 
water-leaving radiance, therefore it is important to apply the NIR correction algorithm (Bailey et al., 2010). 
Unfortunately, even for these carefully selected spectral bands, the 787nm channel must still be corrected 
for contamination by oxygen. Due to the application of spectral smoothing in the calibration processing of 
HICO, the effective spectral bandwidth for wavelengths > 745 nm is 20 nm, so the 787 nm is affected by 
the Oxygen A-band (759-770 nm) (Gao et al., 2012). Oxygen also absorbs light at the Oxygen B-band (688 
nm) and slightly at ~630 nm. Water vapor slightly absorbs the visible light at selective wavelengths and is 
a dominant absorber in the SWIR spectral range, with some notably strong absorption features in the red 
part of the spectrum where we wish to retrieve the water-leaving radiance or its normalized “remote sensing” 
reflectance, R,-s(A): 
f 


Rys(A) = Q™Y) x Lw(d)/Ea(a) (2) 


where E4(A) is the downwelling solar irradiance calculated as: 


Ea(X) = Tyo(A) x to(A) x Fo(A) x fo(A) x cos8o (3) 


Tgo(A) and to(A) are the gaseous transmittances and transmittances of solar irradiance from Sun to surface, 
respectively. Fo(A) is the mean extraterrestrial solar irradiance, f(A) is the Earth-Sun distance correction 
for the time of the observation, and 69 is the solar zenith angle. In Eq. (2), &(A) is the bidirectional 
reflectance correction factor to account for the inhomogeneity of scattered light below the ocean surface and 


the effects of transmittance through the air-sea interface (Morel et al., 2002). 


It is particularly important to retrieve accurate R, (A) across the visible spectral range when applying 


results to algorithms that utilize the spectral shape or spectral derivatives to derive bio-optical properties 
such as phytoplankton functional types (“IOCCG (2014). Phytoplankton Functional Types from Space,” 
2014; Lubac et al., 2008; Ryan et al., 2014). Residual spectral features in R,,(A) resulting from inaccurate 
compensation for water vapor absorption will degrade our ability to retrieve these advanced bio-optical 


properties. 


2.2. Absorbing gases correction 


The gaseous absorption term, T,(A), in Eq. (1) is calculated from the ATREM program (Gao and Davis, 
1997), which calculates the Line-by-Line (LBL) transmittance of seven gases: water vapor (H2O), carbon 
dioxide (COz), ozone (O3), nitrous oxide (N2O), carbon monoxide (CO), methane (CH,), and oxygen (O2) 
for the spectral region 400-3000 nm. The ATREM algorithm was integrated into L2GEN to provide hyper- 
spectral compensation of the water vapor absorption for the atmospheric correction process. Concentrations 
of other gases such as Oz, O3, and NOz, were determined using ancillary sources (Dee et al., 2011; Derber et 
al., 1991; Kanamitsu et al., 2002). The ancillary information is ingested into L2GEN to calculate the trans- 
mittance along the two-way path from the Sun to the surface to the sensor. In L2GEN, ATREM calculates 
the transmittance using the HITRAN 2012 database at 0.05 cm! wavenumber spectral resolution, which is 
down-sampled to a sensor spectral resolution. ATREM has the capability to estimate the column water vapor 
amount (CWV) (given the water vapor volume mixing ratio (VMR) profile) and correct for the water vapor 
absorption along the radiant path. It does so by using a 3-band ratio technique utilizing two atmospheric 
window channels around one absorbing water vapor band. The trough in the TOA reflectance at the water 
vapor band relative to the two window bands is correlated to the water vapor amount along the path that 
the light traveled. The 3-band ratio algorithm assumes that the surface reflectance is spectrally monotonic 
and that any surface spectral features within the spectral windows can lead to erroneous correction for water 
vapor transmittance losses. Fortunately, ocean reflectance is generally monotonic in the red/NIR part of 
the spectrum, except in extreme bloom conditions and floating vegetation (i.e. Cyanobacteria blooms where 
there is a peak in reflectance around 710 nm) (Kudela et al., 2015). 

For HICO, the spectral range contains three prominent absorbing water vapor bands, 725, 825, and 940 
nm. Due to the radiometric calibration issues, the 825 and 940 nm bands are unavailable for determination of 
water vapor correction with HICO, thus the 725-nm band is used to estimate the water vapor transmittance 
along the path. The two window channels are selected to be 705 and 740 nm since they are minimally 
impacted by other gases absorption (i.e., oxygen), and spectrally close, thus reducing the uncertainties 
associated with a non-monotonic surface reflectance (Gao et al., 1993a). The water vapor amount estimated 


from these 3-bands is used to calculate the spectral transmittance along the path, which is combined with 


the derived transmittance effects of the other gases to form T, (A) of Eq. (1). A detailed discussion about 
band selection is in section 2.4. The transmittance of NOz, O3, HzO, O2 absorbing gases, and the total 


transmittance, modulated by HICO’s spectral response function, are shown in Fig. 1. 
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Figure 1: The gases transmittance for a tropical atmosphere for HICO channels with ozone concentration of 277 DU and 
column water vapor of 3.3 cm. 


The LBL calculations for every pixel in the scene (i.e. varying Sun and viewing geometry) are compu- 
tationally intensive. A number of methods have been developed to improve the computational efficiency of 
the LBL approach (Gao and Kaufman, 2003). In this work, we adapted the k-distribution method, which 
involves the grouping by strength of the gas absorption coefficients within a frequency or wavelength interval 
(Kato et al., 1999; Lacis and Oinas, 1991). The interval is selected either based on the sensor spectral 
bandwidth or some set of subintervals with far fewer points to represent the spectral absorption relative to 
the LBL computations. A carefully selected binning or interpolation technique is needed to accurately rep- 
resent the k-bins. Since ranking the LBL absorption coefficients in every subinterval will provide a smooth 
cumulative frequency distribution function (CDF), a 7-point interpolation scheme allowed us to model the 
coefficient accurately, thereby reducing thousands of points in every interval to only 7. The transmittance 


for every interval is calculated as follows: 


and 


Kg) = | BS de (5) 


The k-binned absorption coefficients are a function of the bins g, where K(g) is calculated from the CDF 


of the LBL absorption coefficient, a(B;), for every band or interval B;. p is the integrated absorbing gas 
density along the path in the atmosphere. The transmittance per band, T(B;) is then calculated using the 
Beer-Lambert law for each of the 19 layers of the atmosphere, and the multiplicity property is applied. To 
demonstrate this mathematical transformation, we show in Fig. 2 the LBL absorption coefficients within 
the 725nm water vapor band with 5.7 nm window, the smooth k-distributed absorption coefficients, and the 


interpolated bins. 
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Figure 2: The LBL water vapor transmittance and the k-distribution coefficients at 725 nm. 


Based on this approach, the water vapor transmittance can be calculated 3 orders of magnitude faster 
than the LBL calculations. Also, this approach gives the flexibility to have dual capabilities of the fast k- 
distribution and slow LBL approach to possibly assess spectral shift calibration using absorbing-gas features 
(Gao et al., 2012). The k-distribution approach was implemented into L2GEN so that changes in the spectral 
response function or the band center of detectors will be adjusted instantaneously in the program without 
the need to recalculate the binning. Although the k-distribution method is an exact solution to dramatically 
reduce the computation, the approximated 7-point interpolation and the limits of subintervals per band can 
lead to some errors. Based on our analysis, the error in the water vapor transmittance can reach 1.5% at 725 
nm in the worst-case scenario, when the amount of the column water vapor is higher than 4 cm. Typically, 
in open ocean conditions, the column water vapor is between 0.5 and 3 cm, thus the error in water vapor 
transmittance is less than 1% (Gao and Kaufman, 2003). In addition, since the spectral sensitivity of HICO 
bands are assumed to be Gaussian distributed (consistent with the calibration algorithm), the tails of the 
distribution will not be represented properly in the k-bin. These issues are sensor specific as long as the band 
filter distribution remains the same. By performing LBL calculations for a small subset of a given scene 
(typically 1 pixel), normalization coefficients are derived to remove any marginal errors introduced from the 


k-distribution method. In that way, the computational time is small, while the error relative to the full LBL 


calculations is reduced to 0.1%. 


2.8. Absorption-scattering coupling effects in the atmosphere 


The transmittance of diffused light, in Eq. (1) and (2), through the atmosphere does not account for 
the coupling between the gases absorption and aerosol + Rayleigh scattering. With increasing scattering 
in the atmosphere, the multiple-scattered radiance is increasingly affected by the absorbing gases along the 
path. Therefore, the treatment of gaseous absorption as predominantly governed by Beer-Lambert law is 
not accurate, especially when the atmospheric signal is dominated by highly-absorbing gases features such 
as the Oxygen A-band and strong water vapor absorption regions. Since we are ignoring these coupling 
effects, it is important to assess the uncertainty associated with the atmospheric correction process. Fig. 
3 shows the TOA reflectance, L;/F'o, for an atmosphere-ocean coupled system for two cases: one includes 
water vapor absorption in the atmosphere, and the other does not, where [4 and F'o are the upwelling TOA 
radiance at nadir and the solar extraterrestrial irradiance, respectively. The solar zenith angle of 10° is used. 
When including aerosol scattering, an exact and efficient vector radiative transfer model is employed for 
the sensitivity study to simulate the TOA reflectance between 600 nm and 800 nm, which is based on the 
successive order of scattering method (Zhai et al., 2009, 2010, 2015, 2017). The atmosphere is assumed to 
be a mixture of molecules and aerosols. The molecular scattering properties are determined by the Rayleigh 
scattering matrix and the optical depth is calculated from the method in Tomasi et al., 2005 (Tomasi et 
al., 2005). The aerosol scattering matrix is determined by the aerosol model developed for ocean color 
atmospheric correction, with 50% of fine mode volume fraction and 80% of relative humidity (Ahmad et 
al., 2010). The aerosol optical depth at 550 nm is 0.2. The aerosol vertical distribution is based on the 
average dust case in Braslau et al., 1973 (Braslau and Dave, 1973). The ocean surface is a rough dielectric 
surface with surface slope distribution determined by the parameterization of Cox and Munk (1954) with 
the wind speed of 6 m/s (Cox and Munk, 1954). There is no ocean water body in this case. For the case 
with water vapor absorption, the water vapor absorption coefficient is obtained from the HITRAN 2012 
database (Rothman et al., 2013). The water vapor volume mixing ratio has a realistic vertical profile that 
is 2.5% at the surface, decreasing to 1% at 2.5 km and almost zero at 10 km, so that the top layers of 
the atmosphere do not have any water vapor absorption. For the case without water vapor absorption, the 
atmosphere is assumed to be conservative with the same scattering properties as the case with water vapor 
absorption. In Fig. 3 (a), the simulations were performed with fine spectral resolution and then averaged to 
an assumed ~ 2 nm square band. The line of “Reconstructed” is calculated by using L;/F for the case with 
no water vapor absorption, multiplied the two-way gas transmittance due to water vapor in the atmosphere, 


which estimates the TOA radiance field by ignoring absorption and scattering interaction. Fig. 3 (b) is the 


10 


percentage difference between the case with water vapor absorption and the line of “Reconstructed”, which 


is a representation of the absorption and scattering decoupling error. 
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Figure 3: (a) is the TOA remote sensing reflectance, Lt /F'o, from 600 to 800 nm simulated for three scenarios: no water vapor 
absorption in green dashed line, in the presence of water vapor absorption coupled with atmospheric scattering in the RT 
simulations in black solid line, and the reconstructed TOA reflectance calculated from the no water vapor absorption 
reflectance multiplied by the two-way gas transmittance due to water vapor in the atmosphere in red dashed line. (b) is the 
absorption-scattering decoupling percentage error. 


Fig. 3 (a) and (b) shows that the decoupling error is large for spectral regions with large water vapor 
absorption. The decoupling error is almost everywhere negative, which means that the “Reconstructed” 
line is smaller than the exact radiative transfer simulation with water vapor absorption included. This is 
because the water vapor profile is concentrated at the bottom of the atmosphere. Many photons in solar 
light are scattered back to space at the top layers of the atmosphere, which does not possess any water 
vapor absorption. These photons do not penetrate to the bottom of the atmosphere, and thus do not 
suffer the water vapor absorption at the bottom of the atmosphere. In the absorption-scattering decoupling 
approximation, no discrimination is used for photons scattered at different layers of the atmosphere, which 
is why the decoupling approach underestimates the TOA reflectance. This effect is stronger for spectral 
regions with stronger water vapor absorption. Fig. 3 (b) shows that this error ranges from 0% to 10% 
for the case studied, which has 2 nm spectral resolution. This decoupling error is highly dependent on the 
spectral resolution, with larger error expected for finer spectral resolution. Decoupling errors can also be 
significant for other remote sensing applications. For example, the aerosol scale height can be significantly 
underestimated over dark surfaces when utilizing the Oxygen A-band (Dubuisson et al., 2001). To augment 
the problem, Dubuisson et al., 2001 suggest to derive a correction factor, based on extensive RT simulations 
for varying aerosol optical thickness and geometries, reducing the algorithm uncertainties by an order of 
magnitude (Dubuisson et al., 2001). The importance of a correction method was first reported in Frouin et 


al., 1990, to improve the estimation of CWV over dark ocean surfaces (Frouin et al., 1990). This method was 


11 


then formulated in Dubuisson et al, 2004, where a correction factor was derived as a function of geometry, 
aerosol optical thickness, and aerosol scale height that reduces the retrieval uncertainties of CWV (Dubuisson 
et al., 2004). For ocean color, decoupling correction is also important for the detection of phytoplankton 
fluorescence near 685 nm, where both the Oxygen B-band and water vapor absorb. Decoupling correction 


methods will be developed and implemented in a future SeaDAS (L2GEN) release. 


2.4. Selection of water vapor bands for HICO 


While our primary goal is to retrieve the gas transmittance for atmospheric correction, with additional 
assumptions we can also retrieve the total column water vapor that is of interest for atmospheric science. 
Using the modified L2GEN, the column water vapor in the atmosphere can be retrieved for a given water 
vapor volume mixing ratio (VMR) profile within the atmospheric column. The sensitivity of unknown column 
water vapor VMR profile, surface reflectance, or the presence of thin clouds can introduce more than 10% 
error in the retrieval of column water vapor (Gao et al., 1993a; Gao and Kaufman, 2003; Kaufman and Gao, 
1992). The uncertainties in removing the spectral absorption features from the TOA radiance, however, 
is less significant (Gao et al., 2000). To understand the algorithm errors related to HICO’s band selection 
for the retrieval and correction of water vapor, we utilized vector radiative transfer simulation program, 
MYSTIC, which is a Monte Carlo based RT code, that is computationally efficient during the absence of 
scattering aerosols (i.e. Rayleigh only atmosphere), available through the open source RT library called 
Libradtran (Emde et al., 2016; Mayer, 2009). Furthermore, since the MYSTIC code is an open source and 
well-benchmarked throughout the atmospheric community, we were attracted towards utilizing MYSTIC 
over the SOS code for the case of high spectral resolution and large number of water vapor simulation 
cases. Meanwhile, a comprehensive benchmark study of the SOS code, including an efficient water vapor 
calculations is currently in progress and has not been published. On the other hand, the SOS case study in 
section 2.3, which is more efficient in handling scattering aerosols, has been compared to the MYSTIC code, 
yielding similar results. The MYSTIC code was used to simulate the polarized TOA radiance for a simple 
Rayleigh atmosphere and absorbing flat ocean, at a very high wavenumber spectral resolution of 1 em7! 
(~0.01nm in VIS). Water vapor is assumed to be the only absorber in the atmosphere, which is coupled 
with the scattering in the VRT simulations. The simulation runs were calculated for 0.3, 0.5, 0.7, 1, 1.3, 
1.5, 1.7, 2, 2.3, 2.5, 2.7, and 3 cm column water vapor at geometries permuted from 10° to 50° solar and 
viewing zenith angle with 10° step, while the relative azimuth was fixed to be 90° composing a total of 300 
cases. The water vapor profile was assumed to be the US standard 1976 and the water vapor absorption 
coefficients were obtained from the HITRAN 2012 database (Anderson et al., 1986; Rothman et al., 2013) 


. The high spectral resolution TOA radiance was then convolved with HICO’s spectral response function. 
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Fig. 4 (a) shows, as an example, the simulated TOA reflectance at 10° solar angle, 30° sensor angle, and 


90° relative azimuth angle. 
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Figure 4: (a) Top of the atmosphere reflectance at 10° solar angle, 30° sensor angle, and 90° relative azimuth angle for 12 
different column water vapor values: 0.3, 0.5, 0.7, 1, 1.3, 1.5, 1.7, 2, 2.3, 2.5, 2.7, and 3 cm simulated for HICO. (b) The 
corresponding water vapor transmittance along the radiance path for 12 different column water vapor values. (c) is an inset of 
the transmittance for wavelengths range from 700 to 1000 nm. 


The TOA reflectance has features imposed by the water vapor absorption along the radiance path. The 
strongest water vapor features are in the 725, 825, and 940 nm regions as shown in Fig. 4 (b). To compensate 
for these features, the water vapor transmittance is calculated from L2GEN and the TOA reflectance is 


corrected by dividing it by the gas transmittance along the path. 
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Figure 5: Scatter plot between the assumed CWV in the VRT simulations and the retrieved CWV using L2GEN using a pair 
of water vapor channels at 820 and 940 nm (red circles), 720 and 820 nm (blue circles), and 720 nm only (black circle). The 
error bar is the standard deviation due to changes in solar and viewing geometries. 


Fig. 5 shows the scatter plot between the assumed CWV in the VRT simulations and the retrieval of 


CWV by L2GEN. In this analysis, we used three modes of retrieval. The first is using the average CWV 
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retrieval using the strongly absorbing 940 nm band and the less sensitive 820 nm band shown in blue circles, 
with the error bar of 1 standard deviation due to changes in the simulated geometries. With this combination, 
the CWV retrieval shows a strong underestimation at higher water vapor concentrations. This is due to the 
combination of a weak ocean signal (in that case ocean surface reflectance) and a strongly absorbing 940-nm 
water vapor feature that leads to loss of sensitivity to further increases in CWV. The retrieval for CWV 
less than 1 cm, however, shows good performance suggesting that the 940-nm channel can be well utilized 
to detect and retrieve small amounts of water vapor in the atmosphere, while its sensitivity saturates for 
larger than 1 cm CWV. The second mode of CWV retrieval shown in green circles, using an average of CWV 
retrievals at both 720 and 820 nm, shows very good performance along the whole dynamic range of CWV. 
Although the retrieval shows a slight bias at small CWV values, the impact of erroneous (biased) CWV at 
low values is less significant on the R,, of the ocean, especially at weakly absorbing bands in the visible range 
of the spectrum. The third mode of CWV retrievals using 720 nm only, shows also a good retrieval along 
the whole dynamic range. Retrievals at low CWV values show less bias compared to the combination of 720 
and 820 nm, while there is a stronger bias at high CWV values. The absolute average percent error (AAPD) 
is 19%, 8.5%, and 9% for retrievals using 820 and 940 nm, 720 and 820 nm, and 720 nm only, respectively. 
In Gao and Kaufman, 2003, their CWV retrieval using MODIS showed a systematic bias relative to both 
a ground-based sunphotometers (AERONET) observations and a smaller bias to a ground-based, upward 
looking microwave radiometer (Gao and Kaufman 2003). In the latter case, the error in CWV retrieval was 
less than 10% which corroborates with the analysis shown here, except in the 820 and 940 nm combination, 
which is not ideal for a wide dynamic range of retrievals. Based on the analysis presented here, it is therefore 
recommended to use either 720 or 820 nm or both for the atmospheric correction of water vapor over oceans. 
For HICO, the 720-nm band is used alone, due to the increased uncertainty for spectral bands above 800 


nm. 


3. Radiometric corrections for HICO 


In addition to the known calibration problems with HICO such as second-order optically diffracted stray 
light effects from the grating system in the NIR-SWIR, smile effects, thermal instability, electronic smear 
and low SNR values, spectral band shifts can also lead to poor performance for the correction of gases with 
strong spectral features, such as oxygen and water vapor. To further reduce the impact of sensor calibration 


errors, a vicarious calibration was developed and applied to HICO. 
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8.1. Center Wavelength Shift 


Since the post-launch spectral calibration of HICO is not well tracked, due to the lack of an onboard 
calibration system, and in order to understand the impact of spectral miscalibration in bands center, an error 
analysis was conducted to evaluate the possible error induced in the gas absorption correction algorithm. The 
spectral miscalibration will predominantly impact the correction for gases with sharp spectral features, such 
as oxygen and water vapor. Fig. 6 shows the error in water vapor transmittance at three prominent water 
vapor features in the VIS/NIR range at 590, 720, and 820 nm. This theoretical analysis assumes 5 nm square 
bands, and it shows the error in the water vapor transmittance of LBL calculations and the un-interpolated 
k-distribution calculations (where they have <0.1% error in transmittance without any miscalibration) with 
shifted from the band centers. It is shown in the figure that the error can be as high as 20% for 2 nm 


band center shift at 720 nm. The error is also a function of the column water vapor, where it is positively 


correlated. 
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Figure 6: The error in water vapor transmittance at 590, 720, and 820 nm as a function of the column water vapor induced by 
a spectral band center shift. 


This analysis shows the importance of evaluating band center shifts and the error associated with it for 
gas corrections. For example, Gao et al., 2012 showed that HICO suffered a band center shift immediately 
post-launch (Gao et al., 2012). The initial spectral shift was -1.72 nm that stabilized to -0.9 nm shift. 
In the analysis by Gao et al., 2012, a spectral matching technique was used to evaluate the spectral shift 
based on the known oxygen feature at 765 nm. The spectral matching technique of the transmittance was 
calculated using the LBL absorption coefficients which are necessary for evaluating sharp gases absorption 
features, as opposed to the pre-computed k-distribution coefficients for a predetermined set of band centers 
and bandwidths. Therefore, the k-distribution is required to be re-computed with the possibly shifted bands 


in order to have a proper gas correction, which is implemented in L2GEN. 
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8.2. Vicarious calibration 

The vicarious calibration of ocean color satellite-borne or airborne radiometers is a system-level calibra- 
tion, wherein a set of spectral gains are derived that adjust the TOA radiances to minimize the differences 
between retrieved water-leaving radiances and high-quality in-situ water-leaving radiance measurements. 
Vicarious calibration is an integral part of ocean color remote sensing because the sensitivity of ocean color 
retrievals to radiometric or algorithmic bias is very high. Although HICO is an experimental hyperspectral 
instrument with known calibration stability issues, a vicarious calibration can significantly improve ocean 
color retrievals. In this paper, vicarious calibration is applied to HICO to minimize residual bias in R,.(A) 
retrievals due to instrument calibration error and/or systematic algorithm error (Franz et al., 2007). All 
NASA-supported ocean color missions have used the Marine Optical BuoY (MOBY) near Lanai, Hawaii, 
which has been continuously operated by NOAA since 1996, as the in-situ calibration source for vicarious 
calibration (Clark et al., 1997). The approach involves two steps: 1) vicariously calibrate the pair of bands 
in the NIR (i.e. 747 and 787 nm) used to estimate the aerosol contribution at the TOA and 2) vicariously 
calibrate the remaining bands at the MOBY site, using the atmospheric correction as calibrated in step 
1. To calibrate HICO’s NIR bands, the only 10 available scenes over the South Pacific Gyre (SPG) were 
used and the vicarious calibration was applied assuming zero water-leaving radiance in the NIR. The clos- 
est atmospheric aerosol models to the 80-aerosol model suite, shown in table 2, estimated from the aerosols 
spectral slope from 443-869 nm (average Angstrém) and relative humidity (RH), were taken from co-located, 
semi-co-incident MODISA aerosol model retrievals, after screening for clouds. When MODISA scenes are 
cloudy, the aerosol models for HICO scenes (Model Angstrém) were assumed based on the ancillary relative 
humidity (RH) and climatological aerosol models of mostly coarse sea-salt particles as in Franz et al., 2007, 
shown in table 2. 

The vicarious calibration of the NIR bands is challenged by low SNR. To improve the SNR, we performed 
an optimization analysis by spatially averaging the ~100 m pixel to larger pixel sizes. The coefficient of 
variation (CV = o/j) for the Angstrém coefficient, which is an indicator of aerosol type, was calculated 
for various binned pixel sizes (Fig. 7). o is the standard deviation and pu is the average of the Angstrém 
coefficient of the entire image.Typically, over the SPG regions, aerosols are of oceanic type (i.e. sea salt 
aerosols), and since these remote regions are not impacted by anthropogenic effects, a spatial homogeneity 
in aerosol type is expected. In Fig. 7, the CV becomes smaller and reaches near-asymptotic beyond 1 km 
pixel size. Based on this analysis, a 1-km pixel resolution was selected for the vicarious calibration analysis. 

Using the vicarious calibration process with the assumed aerosol models shown in Table 2, a pair of gains 
was derived with 1 km pixel size. The derived gain factors are 1.038 and 1 at 747 and 787 nm, respectively. 


To vicariously calibrate the visible and other NIR bands, nearly 100 HICO scenes were collected over the 
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Table 2: Selected SPG scenes for NIR vicarious calibration 


HICO Date Time Co-located MOIDSA MODISA _ Average Model 

file name (UTC) file name average RH(%) Angstrém 

Angstrém applied to 

HICO 

H2010079203605  2010- 02:04 A2010032014000 & cloudy 72.91 0.6848 

02-01 A2010032014500 

H2010122031458 — 2010- 08:31 A2010063010000 cloudy 76.67 0.7651 
02-19 

H2010157020639 = 2010- 12:21 A2010061011000 0.9028 78.68 0.8527 
03-03 

H2010162222953  2010- 12:23 A2010061223500 cloudy 77.85 0.7877 
03-02 

H2010216003529 = 2010- 04:03 A2010081235000 0.8264 79.53 0.8033 
03-22 

H2011099205338  2011- 08:05 A2011044013500 cloudy 73.42 0.7651 
02-12 

H2011231021638 2011- 06:22 A2011085013000 cloudy 73.10 0.7651 
03-25 

H2012103031211  2012- 08:31 A2012042011500 cloudy 75.61 0.7651 
02-10 

H2012151211730 2012- 02:12 A2012054014000 cloudy 71.52 0.6848 
02-22 

H2012227014530  2012- 04:14 A2012083011000 cloudy 71.51 0.6848 
03-22 


MOBY site. After filtering for cloudy scenes, and limiting co-incidence of HICO’s overpass with the in-situ 
measurement to +/- 3 hours, 17 scenes were selected to perform the vicarious calibration. Fig. 8 shows the 
derived median and standard deviation of the hyperspectral gains (the full list of gain values is provided 
in the Appendix). Due to the calibration problems in the UV and NIR bands > 800 nm, the gains were 
unstable (due to secondary diffracted light effects on the detector and lack of thermally stable system) and 
deviate dramatically from 1, while the average gains between 400 - 800 nm showed much better consistency, 
with values generally within 5% of unity. We have also made sure to examine the gain coefficients in the 
water vapor absorption spectral bands < 800 nm, and we noticed that there are no or little water vapor 
features in the gain coefficients. We believe that the water vapor correction was done properly even before 


the vicarious calibrations as discussed in our theoretical analysis in section 2.4. 


4. Results and discussion 


The atmospheric correction algorithm and the gains derived from the vicarious calibration process were 
applied to the HICO observations to retrieve hyperspectral R,,. As a verification of system performance, 
the approach was applied to all HICO scenes available over the MOBY site. In Fig. 9, R,, derived from 
HICO after the atmospheric correction process, with and without applying the vicarious gain factors, are 


compared to MOBY’s in-situ R,, optically integrated to HICO’s spectral response function. Also shown is 
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Figure 7: The CV and normalized CV of the Angstrém coefficient for different pixel sizes for 10 HICO scenes collected at the 
SPG. The black vertical dashed line is at the 0.5 km aggregated pixel, while the black vertical solid line is at 1.1 km 
aggregated pixel. 
co-incident R,, from MODISA, when available. It is clear that the R,, match-ups from HICO are improved 
after applying the vicarious calibrations, showing a good agreement with both in-situ MOBY and MODISA 
retrievals. HICO’s R,. also does not contain any features from the absorbing gases (i.e., negative reflectance 
at the 720-nm and 820-nm water vapor bands), including at the water vapor bands, emphasizing that the 

gaseous compensation process is performing well. 

The hyperspectral comparison of HICO and MOBY R,, is very good, except at wavelengths > 800 nm 
where the R,., deviates from the expected zero reflectance possibly due to the second order light contribution 
from the visible to the NIR and thermal instability, effects for which the vicarious calibration can only account 
for mean systematic bias impacts. The average root mean square error (RMSE) for calibrated HICO is 0.001, 
while before calibration it is 0.0024. The improved NIR-band calibration, which that determines the aerosol 
contribution for the atmospheric correction, reduces the bias in the visible spectrum. Not shown here, 
radiometric performance of HICO in the UV wavelengths is also unreliable, but MOBY does not provide 
measurements in the UV to support vicarious calibration in that spectral range. Overall, the MODISA 
shows very good agreement in R,, with MOBY, as expected, since the vicarious calibration was performed 
at the same site. In addition, and unlike HICO, MODISA’s calibration includes continuous monitoring 
and correction for temporal instabilities, using on-orbit calibration measurements of the Sun and the Moon 


(Franz et al., 2007; Meister et al., 2010; Sun et al., 2007; Xiong et al., 2010). 


4.1. Validation 


To provide a comprehensive validation of the atmospheric correction and vicarious calibration, we per- 


formed a matchup analysis between HICO R,., retrievals and NASA’s SeaBASS and AERONET-OC global 
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Figure 8: The hyperspectral median and standard deviation of gains derived from the vicarious calibration of 17 HICO scenes 
at MOBY. 


in-situ dataset (Werdell et al., 2002; Zibordi et al., 2009). Hyperspectral match-ups were not possible because 
no successful match-ups were found within a +/- 4-hour window of the data available, therefore we had to 
rely on creating match-ups with the more abundant multi-spectral data. Figure 10 shows the map of in-situ 
data points around the world. Most of the in-situ observation are located at AERONET-OC sites, since it’s 
more feasible to target a fixed platform than a moving boat, especially when HICO acquired limited number 


of scenes and therefore, there are only very few SeaBASS data points. 


The R,, match-ups were performed for 6 bands in the visible spectral range near 412, 443, 490, 547, 
551, and 667 nm. Fig. 11 shows the scatter plot comparison between HICO’s R,., retrievals and in-situ R,-s 
within a +/- 2-hour window of the overpass. The figure shows the comparison with and without applying 
the vicarious gain factors. As expected, applying the gain factors shows an improvement in the match-ups, 
especially in the blue-green wavelengths, where they are more sensitive to atmospheric correction. The in- 
situ data points represented measurements from various water conditions, from few clear water to mostly 
coastal regions, indicating that the vicarious calibration coefficients can be applied to globally acquired 
scenes. Table 3 shows a summary of the statistical indicators for the match-ups. Additionally, we performed 
a nonparametric hypothesis test using the Mann-Whitney U-test in order to evaluate the impact of the 
vicarious gains on the validation dataset. The p-values shown in table 3 indicate the significance of the 
vicarious gains for each band. The low p-values indicates the rejection of the null hypothesis (i.e. without 


vicarious gains), where the common rejection value is p<0.05. 
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Figure 9: The hyperspectral R;; matchups between in-situ MOBY and HICO retrieval with and without applying the vicarious gain factors. The Figures title shows the scene 
ID from the NASA’s oceancolor webpage. 
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Figure 10: Location of global in-situ SeaBASS and AERONET-OC validation match-ups with HICO shown in red circles. 
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Figure 11: The scatter plot showing the multispectral R-s; comparison between global in-situ SeaBASS and AERONET-OC 
dataset and HICO retrieval. 


4.2. Case studies 


HICO’s capability to perform ocean color retrieval has been improved with hyperspectral atmospheric 
correction and vicarious calibration. The study case comparisons in this section between HICO and MODISA 
show the potential of HICO’s ocean color applications spatially aggregated of TOA at 1-km pixels, but not 
down-sampled. Two locations were selected for comparison: 1) the Chesapeake Bay area, east coast of 
the United States; a highly complex and productive estuary with large anthropogenic influences on both the 
ocean and the atmosphere, and 2) the Gulf of Mexico coastal region off the northwest coast of Florida, United 
States, where phytoplankton blooms are frequent. In both cases, the atmospheric correction is challenging, 
as non-negligible ocean reflectance in the NIR can be significant. Fig. 12 (a) and (b) are the RGB true color 
image composites of the TOA radiance at 461, 553, and 639 nm for the Chesapeake Bay scene, and Gulf of 


Mexico scene, respectively. 
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Table 3: A summary of statistical indicators to in-situ match-ups of HICO with in-situ R;s with and without vicarious 
calibration. The values in parenthesis are without vicarious calibration. 


Wavelength (nm) 412 443, 490 547 551 667 
N 33 45 AT 45 46 42 
R? 0.71 0.73 0.82 0.77 0.77 0.74 

(0.64) (0.67) (0.75) (0.72) (0.72) (0.74) 

RMSE 0.0038 0.0027 0.0054 0.0011 0.0013 0.0004 
(0.019) (0.0112) (0.0104) (0.0091) (0.0084) (0.0023) 

slope 0.95 1.20 1.06 1.03 1.03 1.14 

(0.98) (1.32) (1.15) (1.08) (1.08) (1.21) 

bias -0.0004 -0.0014. —--0.0013.—--0.0004_~—--0.0004_ ~—-0.0003 
(0.0034)  (-0.0001) (0.0005) (0.0008) (0.0007) (0.00001) 

p-value 9x10- 1xi10-7 9x10-> 0.002 0.0032 0.0014 
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Figure 12: (a) true color image composite of HICO’s TOA radiance at 461, 553, and 639 nm of the Chesapeake Bay. (b) true 
color image composite of HICO at the Gulf of Mexico. 


The results in Fig. 13 show the map and histogram comparisons of [Chl] and R,, at 443, 488, and 560 
nm retrieved at the Chesapeake Bay by the NASA’s standard three-band empirical band ratio algorithm 
(OC3) (O'Reilly et al., 1998, 2000) as compared to co-incident, co-located MODISA OC3 retrievals. The 
time difference of the acquisition between the two scenes is nearly 2.5 hours. 

Notably, the average agreement in [Chl] between the two scenes is very good, with HICO showing more 
spatial details due to the higher native resolution as compared to the 1-km MODISA pixel. The agreement in 
[Chl] of the highly dynamic inland estuaries is less noticeable, which can be due to the time difference between 
MODISA and HICO scene acquisition. Spatial anomalies in HICO’s [Chl] image can be apparent such as the 
feature of increased [Chl] across the bay area. That anomaly is attributed to erroneous [Chl] retrieval due to 
the increased contribution of reflectance backscattered from the Lucuius J. Kellam Jr. Bridge-Tunnel across 
the mouth of the Bay. The high spatial resolution in HICO was able to discern the man-made structure, 
whereas MODISA’s image shows much less sensitivity due to the high spatial averaging over 1 km pixels. 


Similar to the Chesapeake Bay analysis, we show the retrievals and comparisons between a HICO scene at 
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Figure 13: A comparison of [Chl] retrieval between MODISA (right map) and HICO (left map) in the Chesapeake Bay and 
histogram comparisons of [Chl] and R,s at 443, 488, and 560 nm, used in the OC3 algorithm. 


a coastal region in the Gulf of Mexico and MODISA. Fig. 14 shows the chlorophyll concentration retrievals 
of HICO and MODISA for the same co-located scene and histogram comparisons of the [Chl] and R,., at 443, 
490, and 560 nm used for the OC3 [Chl] algorithm. The time difference between the two scenes is nearly 1.5 
hours. The spatial variability and magnitude of the [Chl] are quite similar for both sensors, and HICO shows 
more details regarding the spatial structure of phytoplankton blooms in the area as compared to MODISA, 
since it’s not down-sampled. The high native spatial resolution advantage in HICO allows a better distinction 
to the blooms’ spatial structure even when HICO was spatially averaged to a 1 km resolution. Across the 
shoreline, the MODISA retrievals show higher [Chl] than HICO at some pixels. This could be due to the 
spatially averaged bright land pixels (i.e., ephemeral impact) and the impact of the sandy benthos of the 
shallow waters. Statistically, there is a strong agreement in both R,, and [Chl]. 

For a more detailed spectral analysis, Fig. 15 (a) shows the hyperspectral R,, retrieved from HICO and 
the multi-spectral MODISA retrievals at three stations in the image from the Chesapeake Bay, while Fig. 
15 (b) is the same as (a) but for the Gulf of Mexico scene. We selected these locations since they represent 
three different water types and [Chl]. In Fig. 15 (a), ST1 is located in the York River, a highly turbid, 
highly productive region. ST2 is located at the mouth of the Chesapeake Bay and ST3 is just outside the 


bay in the Atlantic Ocean. In Fig. 15 (b) ST1 is located in eastern Choctawhatchee Bay, ST2 is in near the 
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Figure 14: A comparison of [Chl] retrieval between MODISA (right map) and HICO (left map) in the Gulf of Mexico and 
histogram comparisons of [Chl] and R,s at 443, 488, and 560 nm, used in the OC3 algorithm. 


coast in the middle of a phytoplankton bloom, and ST3 is further away from the coast in a relatively clearer 


Gulf water. 


The agreement between MODISA and HICO R,., retrieval is very good for the three locations, indicating 
good consistency in algorithm performance regardless of the water type. Fig. 15 also demonstrates the 
spectral features that HICO can resolve as compared to the multi-spectral MODISA. In Fig. 15 (b), the 
agreement between MODISA and HICO is very good as well, similar to the Chesapeake Bay case. The 
retrieval at ST1 shows unphysically negative R,, values in the blue channels for HICO and for MODISA at 
412 nm. The unphysical value may be attributable to the problematic atmospheric correction over highly 
turbid waters that nullify the black pixel assumption in the NIR and therefore overestimate the aerosol 
contribution. Although, the NIR correction algorithm of non-black pixel was applied to the image, the 
bio-optical mode in the algorithm doesn’t work in those extreme turbid conditions (Bailey et al., 2010). 
It is also possible that the scene is contaminated with absorbing aerosols, which cannot be detected with 
HICO or MODISA spectral capabilities, or impacted by adjacency effects or benthic effects. Such cases of 
high turbidity or complexity are extreme and represent less than a few percent of the global ocean. These 


scenarios can be potentially addressed with some increasing spectral capabilities by utilizing UV channels, 
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Figure 15: (a) and (b) is the comparison of the spectral R;; between MODISA (dashed lines) and HICO (solid lines) 
retrievals at three locations in the Chesapeake Bay, and the Gulf of Mexico, respectively. 


or with polarimetric observations, as planned for the PACE mission. 


5. Conclusion 


Hyperspectral atmospheric correction was developed and integrated into NASA’s multispectral atmo- 
spheric correction code, L2GEN. The extension of the algorithm allowed the compensation of absorbing 
gases in the atmosphere, mainly for the water vapor, and the correction of Rayleigh and aerosol scatter- 
ing over the visible hyperspectral range. The algorithm was tested and validated with the HICO sensor, 
where retrievals of the hyperspectral R,, and [Chl] were compared to in-situ measurements and MODISA 
retrievals. HICO’s R,, retrievals didn’t show spectral features attributed to the absorbing gases, indicating 
a proper compensation in the atmospheric correction. In order to provide ocean color products, the vicarious 
calibration process was applied to HICO based on MOBY’s in-situ measurements. Hyperspectral vicarious 
gains were derived and applied, and R,, retrieval performance was assessed relative to in situ match-ups 
from SeaBASS and AERONET-OC. We also showed two study cases where HICO scenes at the Chesa- 
peake Bay and the Gulf of Mexico are compared to MODISA retrievals. The qualitative and quantitative 
indicators showed a very good agreement with MODISA and in-situ measurements, providing confidence in 
the hyperspectral R,, and derived ocean color products retrieved from HICO observations. This study also 
provides support for atmospheric correction and ocean color algorithm development for the PACE mission, 
where HICO data can serve as a proxy data for hyperspectral algorithm development and testing. The 
full hyperspectral AC algorithm is currently available through SeaDAS version 7.4 and is publicly available 
through https: //oceancolor.gsfc.nasa.gov/. HICO Level-2 scenes will also be available for the community in 


the near future. 
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Table Al: HICO’s hyperspectral gains derived from the vicarious calibration 


Wavelength Gain Wavelength Gain 


(nm) factors (nm) factors 
404 1.042 605 1.038 
409 1.012 610 1.038 
415 1.009 616 1.038 
421 1.025 622 1.050 
426 1.034 627 1.060 
432 1.017 633 1.063 
438 1.019 639 1.063 
444 1.028 645 1.058 
450 1.035 650 1.057 
456 1.046 656 1.046 
461 1.051 662 1.044 
467 1.044 668 1.043 
473 1.041 673 1.041 
479 1.041 679 1.041 
484 1.023 685 1.023 
490 1.013 690 1.026 
496 1.020 696 1.063 
501 1.013 702 1.076 
507 1.013 708 1.072 
513 1.000 713 1.069 
519 0.991 719 1.058 
524 0.993 725 1.059 
530 0.998 731 1.066 
536 1.006 736 1.059 
542 1.014 742 1.052 
547 1.020 747 1.038 
553 1.021 753 1.022 
559 1.015 759 0.986 
564 1.016 765 0.973 
570 1.028 771 0.995 
576 1.038 776 1.000 
582 1.039 782 1.005 
587 1.036 787 1.000 
593 1.033 794 0.992 
599 1.038 799 0.981 
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